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ABSTRACT 

We investigate stationary gas flows in a fixed, rotating barred potential. The gas 
is assumed to be isothermal with an effective sound speed c s , and the equations of 
motion are solved with smoothed particle hydrodynamics (SPH). Since the thermal 
energy in cloud random motions is negligible compared to the orbital kinetic energy, 
no dependence of the flow on c s is expected. However, this is not the case when shocks 
are involved. 

For low values of c s an open, off-axis shock flow forms that is characteristic for 
potentials with an inner Lindblad resonance (ILR) . Through this shock the gas streams 
inwards from X\ to x 2 -orbits. At high sound speeds the gas arranges itself in a different, 
on-axis shock flow pattern. In this case, there is no gas on a^-orbits, demonstrating 
that the gas can behave as if there were no ILR. The critical effective sound speed 
dividing the two regimes is in the range of values observed in the Milky Way. 

We give a heuristic explanation for this effect. A possible consequence is that star 
formation may change the structure of the flow by which it was initiated. Low-mass 
galaxies should predominantly be in the on-axis regime. 

A brief comparison of our SPH results with those from a grid-based hydrodynamic 
code is also given. 
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1 INTRODUCTION 

Gas kinematic observations in the best-studied barred galax- 
ies imply highly non-circular gas motions. Some examples 
are NGC 5383 (Duval & Athanassoula 1983), NGC 7496 
(Pence & Blackman 1984b), NGC 1365 (Teuben et al. 1986), 
and our Galaxy (Liszt & Burton 1980, Binney et al. 1991). 
These non-circular motions are believed to be associated 
with shocks; large velocity gradients in support of this inter- 
pretation have been observed in NGC 6221 (Pence & Black- 
man 1984a), NGC 1365 (Lindblad & Jorsater 1987) and 
NGC 3095 (Weiner, Williams & Sellwood 1993). The shocks 
can be traced on optical images by the dust lanes that mark 
swept-up dense gas at the leading edges of the bar (Pence & 
Blackman 1984a, Athanassoula 1992). In IC 342, millimeter 
observations of the molecular gas distribution and velocity 
field suggest that the gas in the shock ridges is falling into 
the nucleus (Ishizuki et al. 1990). Some estimates for the 
resulting mass infall rates are ~ 0.01 — O.lMo/yr in our 
Galaxy (Gerhard & Binney 1993) and ~ 4M /yr in NGC 
7479 (Quillen et al. 1995). The infalling gas may settle on an 
inner ring and fuel starburst activity (Ishizuki et al. 1990, 
Kenney et al. 1992). 

Because the effective temperature (cloud velocity dis- 
persion) of the gas is much smaller than any orbital veloci- 



ties, it is generally a good approximation to think of quasi- 
stationary galactic gas flows in terms of the periodic orbits 
in the underlying gravitational potential, so long as these pe- 
riodic orbits do not cross. Many hydrodynamic simulations 
using a variety of numerical techniques have confirmed this 
(Sander & Huntley 1976; Schwarz 1981, 1984; Combes & 
Gerin 1985, Habe & Ikeuchi 1985; van Albada 1985; Mul- 
der & Liem 1986; Athanassoula 1992; Friedli & Benz 1993; 
Jenkins & Binney 1994). When the periodic orbits do cross 
or intersect, then pressure or viscous forces must always be- 
come important and ensure that the gas streamlines match 
together to form a well-defined flow pattern. This typically 
happens near the orbital resonances of the system, and it is 
here that shocks are found to form. 

The geometrical structure, the extent, and the strength 
of these shocks depend on a variety of external parameters. 
Most important is whether the gravitational potential has 
an Inner Lindblad Resonance (ILR) and hence a family of 
X2-orbits , but also the strength and axis ratio of the bar 
are significant parameters (Roberts, van Albada & Huntley 
1979, Sanders & Tubbs 1980, Athanassoula 1992). 

We have begun a study aimed at better understanding 
these gas flows, and using them as tracers for the dynam- 
ical properties of our Galaxy and barred galaxies in gen- 
eral. We have used a two-dimensional smoothed particle 
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hydrodynamics (SPH) method, based on the SPH-code of 
Steinmetz & Midler (1993) which Matthias Steinmetz kindly 
made available. The extent to which the interstellar medium 
(ISM) in galaxies can be modelled by simple gas dynamics 
has been discussed recently by Sellwood & Wilkinson (1993) 
and Binney & Gerhard (1993). In the present study, we have 
omitted the self-gravity of the fluid, treating it as a tracer; 
the restriction to two dimensions is in order to increase the 
resolution. Later, the SPH method will allow us to study the 
effect of the gas' self-gravity, particularly in the high-density 
regions that generally form in these flows. 

In the course of this work, we have found that the shock 
properties in barred potentials also depend on internal pa- 
rameters used to model the gas; specifically, on the sound 
speed if the gas is assumed to be an isothermal fluid. This 
paper reports on these results. We first describe the mass 
model, initial conditions, etc. (Section 2) and give some de- 
tails on the numerical method (Section 3). In Section 4 we 
show that both off-axis and on-axis shocks may form in a 
potential with an ILR, study the influence of some relevant 
parameters, and briefly compare our results with those from 
grid-based simulations. Finally, in Section 5 we discuss some 
implications of our main result. 



2 DESCRIPTION OF THE MODELS 

In this paper we consider gas flows in galaxy models with 
a given gravitational potential for the stellar component, as 
detailed in Section 2.1. Our assumptions concerning the ini- 
tial conditions, the hydrodynamics of the gas, and a simple 
gas recycling law for mimicking star formation are given in 
Sections 2.2-2.4, while the numerical method (Smooth Par- 
ticle Hydrodynamics, SPH) is briefly described in Section 3. 



2.1 Mass models for the stellar component 

In order to be able to compare our results from SPH with 
those from a hydrodynamic grid method, we have adopted 
the analytic family of models from Athanassoula (1992). 
These models include a bulge, disk, and bar component. 
The bulge follows a Hubble profile 



p(r) = p b (l+r 2 /r 2 b )- 3/2 



(1) 



and the stellar disk is a simple Kuzmin-Toomre disk with 
surface density 

RdM d 



2tt 



-3/2 



(2) 



The parameters are chosen such as to give an approximately 
flat rotation curve: R d = 14.1 kpc, M d = 2.3 x lO n M , 
rt = 0.33 kpc, and p b = 23.6 M /pc 3 . 

The bar component is modelled as a simple Ferrers el- 
lipsoid with density 



P(X) = 



Mi- 
o 



for g 2 < 1 
elsewhere, 



(3) 



with g 2 = X 2 /a% + (Y 2 + Z 2 )/b% and parameters hf = 
5 kpc, &f = 2 kpc. We have considered two values for the 
bar mass: one is specified by setting pF = 0.4476 M Q /pc 3 
(Model A, corresponding to Athanassoula's model 001), the 
other has one half that mass (our standard Model S). To 



Figure 1. Circular rotation curves in the equatorial plane: total 
(full and dotted lines) for Models A and S, and contributions from 
the bulge (long-dashed), Ferrer's bar (short-dashed for Model A, 
dotted for Model S), and disk (dash-dotted). 
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0.6 
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25 
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1.2 


0.6 


S-5 


30 


no 


1.2 


0.6 


S-6 


15 


yes 


1.2 


0.6 


S-7 


20 


yes 


1.2 


0.6 


A-l 


10 


no 


1.8 


0.4 


A-2 


15 


yes 


1.2 


0.6 


A-3 


20 


yes 


1.2 


0.6 



Table 1. Parameters for the models computed. For all runs we 
used N ft: 2 10 4 particles in 2D and imposed point symmetry, 
giving a mean resolution lengths of h < 0.1 kpc. All simulations 
are rotating clockwise. Table columns give (1) Model designation, 
(2) effective sound speed of the gas, (3) whether gas recycling 
was used, (4) constant determining the individual particle sizes, 
cf. Section 3, and (5) time since start of the calculation for which 
results are reported. 



place the Ll-Lagrange point at t-li = 6.0 kpc we chose a 
pattern speed of Q p = 35.3 km s _1 kpc -1 for the Model 
A bar and Q p = 32.0 km s" 1 kpc -1 for the Model 5* bar. 
This corresponds to bar rotation periods of 0.174 Gyr and 
0.192 Gyr, respectively. Then from the m = 0-component 
of the potential, the corotation radius is -Rcr = 5.8 kpc in 
both cases. 

The rotation curves for both models and the contri- 
butions of the individual components (from the respective 
m = 0-components of the gravitational force in the equato- 
rial plane) are shown in Fig. [[]. 



2.2 Initial conditions and symmetry 

In the present two-dimensional calculations the gas is ini- 
tially set up on circular orbits with constant surface den- 
sity Eo = 1 Mq pc -2 . Outside some radius r cut beyond 
corotation the density is set to a constant, i.e., the pressure 
gradients are set to zero. 

The non-axisymmetric part of the potential is turned 
on gradually within a time r on equal to approximately one 
half the rotation period of the bar. After continuing the 
calculation for some time t given in Table 1 the gas flow is 
nearly stationary in the rotating frame. 

In these calculations we have assumed point symmetry 
with respect to the origin. This effectively doubles the num- 
ber of particles, and increases the resolution by a factor %/2- 
We have checked that the results in simulations without the 
assumed point symmetry are essentially identical. 



2.3 Gasdynamics 

For the gas we have assumed an isothermal equation of state 
with an effective sound speed c s , representing the random 
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motions in the interstellar cloud medium. This approach is 
consistent with the ISM model of Cowie (1980), who argued 
that the cloud fluid can be treated like an isothermal gas if 
the clouds have an equilibrium mass spectrum. The equlib- 
rium is assumed to be maintained by a steady supply of 
small clouds by supernovae. Using SPH as described in the 
next section we are thus solving the Euler equations 



for effective sound speeds of c s = 10 — 30kms . 



(4) 



2.4 Gas recycling 

Because the gas loses angular momentum to the stellar 
bar by gravitational torques, a substantial accretion rate 
towards the centre develops. In run S-l, for example, at 
a radius of 1 kpc, the accretion rate is approximately 
7.6 x 10 6 M /Gyr, and in S-5, it is ~ 1.9 x 10 7 M /Gyr, ris- 
ing continuously with sound speed between these two cases. 
For comparison, for our initial surface density of IMq/pc 2 
the total gas mass inside corotation is 1.06 x 10 s Mq. Thus 
after ~ 5 — 15 Gyrs all of the gas would have been accreted 
to the centre. To avoid this, we have in some of our models 
introduced a gas recycling law which takes away gas in high 
density regions and adds material uniformly over the disk: 



dE/dt = a£o ~ aS 



(•») 



This may crudely be thought of as simulating star for- 
mation in dense regions and mass loss from preexisting 
bulge and disk stars. The star formation rate a is taken 
to be 0.3Mq 1 pc 2 Gyr _1 . Equation (|^) includes a quadratic 
Schmidt (1959) law and a constant source term, but the 
functional form was mainly chosen to be able to compare 
with the grid-based simulations of Athanassoula (1992). 



3 NUMERICAL METHOD 

We have used a two-dimensional SPH method to solve the 
hydrodynamical equations for the gas flow in the galaxy 
equatorial plane. Restricting the calculation to two dimen- 
sions increases greatly the resolution we can achieve, and 
this is important in the present problem for resolving the 
shocks that dominate the evolution. Since we are interested 
in isothermal shocks, we are assuming implicitly that the gas 
can cool rapidly. We then expect that the gas can expand 
vertically with speeds only of the order of the sound speed. 
Since this is much smaller than the dynamical velocities, 
vertical motions cannot affect the shock front even if this 
has finite width (such as one or two cloud mean free paths). 
The gas will also only marginally be able to expand out of 
a three-dimensional layer within a dynamical time; hence 
the expansion of the post-shock gas can not lead it too far 
above the downstream material that it will meet next. Thus 
we expect the two-dimensional description to be a reason- 
able approximation. Notice that it is not possible to test 
this rigourously with three-dimensional simulations of this 
particle number, unless anisotropic particles are used; for to 
resolve the disk vertically, one needs to make it substantially 
thicker than in reality, which prolongs the vertical escape 



times. We have done some such three-dimensional simula- 
tions; in these the observed gas flows were indeed similar to 
the two-dimensional ones. 

As there are many variants of SPH (see the reviews of 
Benz 1990, Monaghan 1992), we give a short overview of 
our code, which is based on the code of Steinmetz&Muller 
(1993). A few details on how the method works will be use- 
ful for interpreting the results reported below. SPH solves 
the equations of motion by a Monte Carlo integration. The 
integration points are not drawn by random, but are the po- 
sitions of the particles, as they happen to be. Ideally the par- 
ticles should be distributed in a glass-like structure to min- 
imize numerical errors. Apart from this the particles have 
no physical meaning. The forces on each parcel of fluid can 
be calculated by smearing out the properties of the particles 
over a few mean particle distances. 

For constant smoothing length h, the smoothing is done 
by folding each field quantity A such as temperature or pres- 
sure with the kernel function W(r) 



(A(r))= / A(r')W(r' -r,h)d 2 x' + 0{h 2 ) 



(6) 



Here A% and Ej are evaluated at the position of particle i. 
Thus the surface density is approximated by 



<£(r)> = ^2m t W(r t ~r,h). 



(7) 



Only structure on scales larger than a few resolution lengths 
h is meaningful. We used the spline kernel 

W(r-r',h) = ^{ 1(2 -uf 

u > 2 



1 < u < 2 



(8) 



with u = \r — r'\/h, and the constant Wo is determined from 
the normalization condition 



W(r,h)d x. 



(9) 



The kernel describes the surface density distribution of one 
particle. 

To increase resolution in high-density regions the par- 
ticles are assigned individual smoothing lengths 

hi = cVrmfSi (10) 

with £ = 1.2 in most cases (see Table 1). In this case, the 
field quantities at particle positions t\ are evaluated with 
the symmetrized Kernel 

Win-r^fu,^) = UW(n-r j ,h i ) + W(r i -r j ,h j )}(ll) 

to ensure momentum conservation. This symmetrized Ker- 
nel is also used in the calculation of the surface density 
Ej = (E(ri)) at particle position n. To evaluate the field 
quantity at an arbitrary position r the original Kernel is 
used. 

We have used the following simple approximation for 
adjusting the particle smoothing length: 



N 



(n) 



(12) 
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Here No = 4n£ 2 is the desired number of neighbours within 
2hi, and A^ n ' the recorded number of neighbours of parti- 
cle i in the n-th time step. This ensures that the number 
of neighbours of a particle is approximately independent of 
time and of its location. In our case with strong shocks the 
original method used by Steinmetz (1993) was found to re- 
sult in a larger number of neighbours in the shock region and 
a loss in resolution. In some cases we also observed large os- 
cillations in the particle size from time step to time step 
when a particle entered the shock. 
The Euler equations in SPH read 



d«i 
~dt 



(13) 



where Qij is an artificial viscosity tensor (Monaghan & Gin- 
gold 1983) 



Qi 



0, 



(n - rj) ■ (vi - Vj) < o 

otherwise, 



(Vi 



(n - rj) 2 + if 



(14) 



(15) 



with a = 1, /3 = 2, r/ = 0.1 hij, hij — (hi + hj)/2, and Ey = 
(Ej + Ej)/2. The quantity [Uj is an approximation for the 
contribution of the particle pair to the local volume 

compression hV-v over one resolution length h. 

The artificial viscosity actually consists of two con- 
tributions. The term with the constant a corresponds to 
the bulk viscosity whereas the term with the constant 
P corresponds to the von Neumann-Richtmyer viscosity. 
Both correspond to a 'viscous pressure' in the Euler equa- 
tion, namely aY,c s \h\7-v\ and /3E(7i V-w) 2 (Benz 1990). The 
von Neumann-Richtmyer viscosity is necessary to simulate 
strong shocks because it guarantees causality by transport- 
ing information with signal speed \f]3\hV-v\ across the 
shock fronts (e.g., the flow is supersonically compressed 
(c s < -ftV'«). The bulk viscosity is used to damp post- 
shock oscillations. 

As reported by Hernquist & Katz (1989), the artificial 
viscosity tensor does not vanish in pure shear flows. There- 
fore, some authors use the modified artificial viscosity tensor 
proposed by Balsara (1995) 



Qij — Qij g (ft + fj) 



with 



fi 



l(V-«)i| 



|(V-«)j| + |(Vxw)i| + O.OOOlCs/fti ■ 



(16) 



(17) 



The effect of this is to reduce the viscosity in regions with 
large values of |(Vxv)i|. However, we have not found this 
viscosity switch useful in the presence of shocks with strong 
shear. In a sheared shock, both the divergence and the rota- 
tion of the velocity field become large, thus the fi decrease 
below unity and the viscosity becomes too small to treat the 
shocks adequately. 

The gas recycling law (jjj]) is realized by adjusting the 



Figure 2. (a) Particle distribution in Model S-2, for time t = 
0.6 Gyr. The outline of the bar is given by the ellipse, (b) The 
velocity field for the same model, (c) Some x\ and X2 orbits useful 
for the interpretation of the gas flow. 



mass of each particle riii at its individual timestep n + 1 
according to 

dm dE 



(71+I) (n) . 1. 

m) — m\ ' + dtj 



dE dt 



= TO C«) +dt .^|j a (E 2 -E 2 ). 
The particle velocities are kept constant. 



(18) 



4 RESULTS 

4.1 Off-axis shocks at low sound speeds 

Fig. ^ shows the distribution of gas particles, the veloc- 
ity field and the underlying orbital structure in Model S-2. 
This model describes an isothermal fluid with a sound speed 
c s = 15kms _ . The velocity field was calculated on a regular 
grid, using the SPH smoothing algorithm, and is shown in 
a rotating frame in which the bar and the gas flow pattern 
are approximately stationary. 

As is well-known, the flow is influenced strongly by the 
families of periodic orbits in the underlying potential and the 
transitions between them. The most important families here 
are the a;i-orbits elongated along the bar and the central 
£2-orbits oriented perpendicular to the bar. The outermost 
orbits shown in Fig. Gr belong to a 4/1-resonant family. See 
Contopoulos & Papayannopoulos (1980) for this notation. 

The velocity field shows that the gas occupies Xi orbits 
in the main bar region and then shifts gradually to X2 orbits 
in the inner parts. In this process it forms a shock along 
the leading edges of the bar; in model S-2 the shock starts 
at w 3.7 kpc and continues inward with a small inclination 
angle with respect to nearby zi-orbits. From a particle point 
of view, the gas particles moving through the shock exchange 
momentum with particles downstream to avoid penetration 
of each other. 

The gas moves inwards alongside the shock and then 
swings around the disk of X2-orbits towards the symmetric 
far-side branch of the shock. The flow follows the a;2-orbits 
some way around the centre, reaches a density maximum 
and opens like a spray as some particles go into the disk, 
whereas others move outward again into the far-side shock. 
With lower orbital energy, they then begin another half- 
revolution. The gas thus loses further orbital energy and 
spirals inwards in a few revolutions until it joins the gas on 
the outer occupied a;2-orbits. The resulting central ring of 
dense gas is often called the a;2-disk. 

Gas coming into the bar region from outside encounters 
a trailing spiral shock wave which occurs at the intersection 
of colliding streams on the outer xi orbits and the 4/1 orbit 
family. In Model S-2 this shock front starts on the major axis 
at w 3.7 kpc where the x\ family meets the 4/1 family. For 
definiteness, we call the two types of shock the 'cusped-orbit 
shock' and the '4/1 spiral shock'. 

Fig. |^ shows the density and the velocities of parti- 
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Figure 3. (a) Surface density, (b) velocity vu s along, and (c) 
velocity v± s perpendicular to 'Slit 1' in Fig. fcj, for all particles 
within one smoothing length h of the slit. This slit intersects the 
cusped orbit shock approximately perpendicularly. The diamonds 
give the surface density £; evaluated at the positions of individual 
particles; the solid line shows the smoothed SPH value along the 
slit computed by integrating over all overlapping particles. Full 
symbols denote particles which are deccelerated in the shock, i.e., 
are compressed supersonically: -feV-u > c s . 



Figure 4. Same as Fig. ^ but for 'Slit 2' through the 4/1-spiral 
shock in Fig. El 



cles within one smoothing length h of 'Slit 1' plotted in 
Fig. [| The left side of Fig. ^ with small s corresponds to 
the low density, pre-shock region inside the cusped orbit. 
The top panel shows the gas surface density through the 
cusped orbit shock along the slit. The points give the sur- 
face density Ej evaluated at the positions of individual par- 
ticles; the solid line shows the smoothed SPH value along 
the slit computed by integrating over all overlapping parti- 
cles [cf. eq. (m)]. As determined from the latter, the density 
jump across the shock is about a factor of two; because of 
resolution effects this is a lower limit. The rise of the density 
from its pre-shock value occurs on a length scale comparable 
to the smoothing length h; here h is approximately 140pc 
upstream, 70pc at the density peak and 90pc downstream. 
The post-shock region (near the peak of the density) ap- 
pears still not resolved; most of the right side of the diagram 
shows densities on streamlines which are quite unrelated to 
the shock (see Fig. ^). 

The lower two panels of Fig. [|show the gas velocities v\\ s 
and v± s parallel and perpendicular to the 'slit' direction in 
Fig. |^a. At the shock surface W|| a is approximately the veloc- 
ity with which a supersonic collision occurs, v± s is approx- 
imately the transverse shear component of the streaming 
velocity. It is seen that v\\ a jumps by ~ lOOkms -1 over the 
same length-scale over which the density jump occurs. The 
full symbols in Fig. [^denote those particle positions at which 
—h V-u > c s — the local criterium for a shock to take place. 
For a simple non-shearing, isothermal shock the Rankine- 
Hugoniot jump conditions predict W|| s , P ro x w ||s, P ost = c s', 
hence for A«ii s — lOOkms -1 and c s = 15kms _1 one finds 
v \\s,pie — 102kms _1 and V|| SlPO st — 2kms _1 for the upstream 
pre-shock and downstream post-shock components, which 
appears approximately correct. The corresponding density 
jump is a factor of 50 and must then be clearly unresolved. 

However, from the lower panel and also from Fig. |^b 
we see that the gas flow in this region has a strong shear 
component and there is likely also a discontinuity in the 
transverse v± s velocity across the shock. Quantitatively, 
A?;x s ~ 180kms _1 across the unresolved shock width. In 
this case, the jump conditions are not as simple (Syer & 
Narayan 1993). 

Fig. ^ shows analogous information for the 4/1-spiral 
shock. This shock is less complicated because shear velocities 
are seen only in the post-shock region, but are essentially 
absent in the shock region itself. The density goes up by a 
factor of four, and the Mach number Vn a /ce also jumps by 
a factor of four. From the Rankine-Hugoniot conditions we 



would expect p pr e/p P ost = 1/M 2 , so the maximum density 
is again unresolved. 

4.2 On-axis shocks at high sound speeds 

So far we have adopted a sound speed of c s = 15 km/s 
which is consistent with observations of the turbulent speed 
of clouds in the ISM. However, the observed values in the 
Galactic disk do vary, so it is interesting to see whether the 
stationary flow pattern depends on c s or not. 

From a theoretical point of view we might expect several 
regimes to occur between two extreme cases: In the limit 
of negligible sound speed every slight compression of the 
gas causes a shock and an isothermal fluid would become 
effectively incompressible. In the other limit of high sound 
speeds no shocks can form at all. 

Indeed, we found that there are two qualitatively dif- 
ferent quasi-stationary gas flow solutions in the potentials 
investigated. The transition between these two configura- 
tions is shown in Fig. [|. At low sound speeds, the gas flow 
is as described in the last section. At higher sound speeds 
(cs — 25 — 30 km s _1 in our standard potential with cir- 
cular velocity « 200kms~ 1 ) the X2-orbits are no longer oc- 
cupied and a cusped orbit shock cannot form; instead, it is 
replaced by a broken shock configuration near the long axis 
of the bar. This structure resembles more the shocks that 
Athanassoula (1992) found to occur in barred potentials 
with no Inner Lindblad Resonance (ILR), while of course 
the potential used still has an ILR. The spiral wave also is 
no longer strong enough to form a 4/1-spiral shock. 

A heuristic explanation for these results is as follows. 
Consider a sequence of models such as in Fig. ^, starting 
with a case with well-developed off-axis shocks. The speed 
with which the spray of gas inside the cusped orbit meets 
the cusped orbit shock decreases outwards (Fig. 2b). Thus 
the shock front should end where the velocity v\\ a in the 
shock frame becomes smaller than the sound speed. When 
the sound speed is increased while keeping the potential con- 
stant, the streaming velocities will to first order remain un- 
changed, so that the shock front must move inwards. This 
means that it must move towards a;i-orbits deeper in the 
potential well and closer to the major axis. The second as- 
pect is the gas inflow along the shock which feeds gas onto 
the X2-ring- As the shock front moves closer to the bar ma- 
jor axis, the inflowing gas will eventually no longer be able 
to settle on X2-orbits. Once a;i-orbits cannot be occupied, 
one of the principal factors instrumental in building up the 
off- axis shock has disappeared. 

From Fig. ||, we conclude that the transition between 
both solutions in terms of the variation of c s is more or 
less continuous. E. Athanassoula (private communication) 
has also found in her grid simulations the existence of 
two distinct gas flow patterns; however, she appears to 
see a very sharp transition at a critical sound speed of 
c s = (16.4 ± 0.1) km s _1 (in the potential model A). It 
is not clear at present precisely which aspects of the numer- 
ical schemes are responsible for this difference; nor whether 
any of these schemes describes the dynamics of the ISM ac- 
curately enough to say how the real ISM behaves in this 
respect. 

An interesting aspect of the flow patterns at high sound 
speed is that the gas inflow still forms a central disk; however 
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Figure 6. An on-axis shock solution switches to an off-axis shock 
flow when the sound speed is changed by hand and the massive 
central xi-disk is removed. Times are specified in Gyr and the 
bar rotation period is 0.174Gyr. 

this now consists of gas on zi-orbits rather than £2-orbits. 
This is only possible because in our Ferrers bar potential the 
x\ orbits are not self-intersecting. If they did self-intersect, 
a still different flow pattern connected with strong accretion 
to the centre would likely result. 

4.3 Influence of initial situation 

We now ask whether the type of shock structure that devel- 
ops in the simulation depends on the initial gas configuration 
as well as on the sound speed. In principle, quasi-stationary 
gas flows with both types of shock structure might exist 
in one and the same gravitational potential, and the sys- 
tem could evolve towards one or the other equilibrium, de- 
pending on its initial state. To test this, we start a simula- 
tion in an extreme initial state where the system is already 
in quasi-equilibrium such that the flow has formed an on- 
axis shock structure, and then see whether it evolves to an 
off-axis shock solution when we reduce the sound speed to 
10 kms -1 . 

When the sound speed is dropped, a shock in the orig- 
inal gas flow must remain a shock. On the other hand, a 
shallow compression wave may suddenly become supersonic, 
and the newly forming shock can change the flow pattern in 
such a way that the original shock structure is modified or 
disappears. When the sound speed is reduced in the equi- 
librium flow shown in the top panel of Fig. [| the main 
shock front immediately moves away from the bar major 
axis in the direction of the position where the main off-axis 
shock would be located in a low-sound speed model. In ad- 
dition, new spiral shocks form near the position of the 4/1- 
resonance, and much of the gas in the inner bar region falls 
towards the centre, building a massive, elongated a;i-disk. 
This xi-disk is massive enough to prevent gas from moving 
onto a;2-orbits, so no off-axis flow can form, and the whole 
configuration swings back towards the bar major axis. The 
resulting quasi-equilibrium is shown in the middle panel of 
Fig. | 

We now turn on gas recycling, which essentially re- 
moves the high density material in the £i-disk. To speed 
up the computation, we use a higher gas recycling constant 
a — 3 Mq 1 pc 2 Gyr -1 than in the other models. This re- 
moves the obstacle in the flow, and after a few more rotations 
the gas now settles on a;2-orbits and an off-axis shock front 
has formed. The structure of the inner £2-disk in the final 
configuration (bottom panel of Fig. ^) is somewhat differ- 
ent from those shown previously because of the large gas 
recycling rate. 

We draw two conclusions from this experiment, (i) If we 
release small quantities of gas in a barred potential, then in- 
dependent of the initial velocity field the flow pattern adjusts 
itself to that which is natural for gas at this sound speed in 
this gravitational potential. For, the velocity field of the gas 
that evolved towards the final off-axis configuration in Fig. ^ 
must have been close to the velocity field typical for on-axis 
flows, (ii) On the other hand, a massive disk set up previ- 



ously is hard to perturb and may prevent the incoming gas 
flow from taking up its natural flow pattern. In real galactic 
disks, such massive disks are expected to form stars and be 
depleted of gas; thus eventually we expect the gas there to 
return to its preferred configuration. 

4.4 Influence of gas recycling and bar mass 

In Fig. Q we compare models with different sound speeds, 
bar masses, and with or without gas recycling as indicated 
in the figure caption and Table 1. 

First we consider the influence of the gas recycling by 
comparing Model S-2 with S-6 and Model S-3 with S-7. In 
both cases the large scale flow configuration does not change. 
However the simulations with gas recycling (S-6 & S-7) look 
more clumpy; we believe this is an artifact of having to intro- 
duce gas particles with different masses in the gas recycling 
algorithm. 

Between the middle panels and the lower panels of Fig. ^ 
the bar mass is doubled and therefore the orbits change as 
well as the resonances. The more massive bar forces more 
elongated xi-orbits (see Fig. ^ below). For the lower sound 
speed the pattern maintains off-axis shocks, but the shock 
front now starts nearer to the end of the bar and is closer 
to the bar major axis. There is also a new large gap in the 
shock front at (x, y) — (2, 2)kpc. 

For the higher sound speed (Model A-3) the flow field 
in the more massive bar potential has changed to an on-axis 
configuration. Thus the critical sound speed which divides 
the two shock patterns appears to decrease for increasing 
bar mass. 



4.5 Comparison of SPH and grid based results 

In Model A-l we have reproduced model 001 of Athanas- 
soula (1992) to test our code against a grid based method; 
in this case, a second order flux splitting code written by 
G. D. van Albada (van Albada & Roberts 1981). Both meth- 
ods solve the Euler equations, thus the same results should 
be expected. But differences may arise due to the statisti- 
cal nature of SPH or the different numerical and artificial 
viscosities used in both codes. Comparing a simulation of 
the same physical problem computed with both numerical 
methods may give valuable clues about their reliability. 

Fig. U shows the particle distribution, gas velocity field, 
and closed orbit structure for Model A-l. The potential is 
the same as used for model 001 of Athanassoula (1992), and 
the initial conditions are set up identically except that the 
disk is truncated at a smaller radius (6kpc as opposed to 
16 kpc in the grid model) . Compared to the simulations dis- 
cussed earlier, we have used a larger smoothing parameter 
(£ = 1.8). The reason for this is that, because of the stronger 
bar, the gas in this potential moves to the center more 
quickly, leaving fewer particles in the shock front in quasi- 
equilibrium. To measure the shock properties displayed in 
Fig. ^, it was then necessary to employ more smoothing. 

The global structure of the flow seen in Fig. ^| is simi- 
lar to that of the grid model (see Athanassoula's Fig. 2 for 
comparison). The point on the major axis where the 4/1- 
spiral shock begins is at the same distance from the centre 
(at ~ 4.3kpc). The cusped orbit shocks are straight in both 
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Figure 5. A sequence of models with increasing sound speed as indicated in the upper left corner of each panel (models S-l to S-5 in 
Table 1). Between 20 and 25 km/s the flow pattern changes as described in the text. The orbital structure in all the simulations is the 
same; some closed orbits are shown in the lower right panel. 



Figure 7. Influence of bar mass, sound speed and gas recycling on the gas flow pattern. Models in the left-hand panels have sound speed 
c s = 15kms — , those in the right-hand panels c s = 20kms~ 1 . Upper two panels show gas flows in the standard potential and without 
gas recycling that mimicks star formation, the middle two panels show gas flows in the standard potential with gas recycling, and the 
lower two panels are for a potential with a bar twice as massive as before and again include gas recycling. See Table I for parameters. 
There is no influence of gas recycling on the flow pattern, but a stronger bar or a higher sound speed can drive the gas flow towards a 
configuration with on-axis shocks. 



Figure 8. Same as Fig. | for Model A-l at t = 0.4 Gyr. This 
model is the same as Model 001 of Athanassoula (1992), and the 
figure should be compared with her Fig. 2. 



Figure 9. Same as Fig. ^ 



for the Slit in 



Fig. | 



simulations, however, in the SPH model they form a smaller 
angle with the bar major axis than in the grid model. 

The shock properties shown in Fig. g are similar to 
Athanassoulas findings (her Fig. 11) when we take into ac- 
count that our and her slits are not at identical positions and 
that the resolution is different. The shape of the function de- 
scribing the gas velocity along the slit is very similar, and 
the maximum upstream value and the downstream veloci- 
ties agree quantitatively (to be able to make the comparison 
with Athanassoula, the velocities in the two lower panels of 
Fig. ^|are measured relative to the inertial frame and include 
a contribution from the bar rotation). The shock is slightly 
less well resolved in the SPH simulation, corresponding to 
a somewhat broader density peak. The peak density is un- 
derestimated in both simulations; in the SPH case it is also 
time-dependent. The value we measure from Fig. ^ appears 
to be higher than that in Athanassoula's Fig. 11 (the units 
in both diagrams are Eo = lM Q /pc 2 ). 

The large scatter in the density plot in Fig. ^ is due to 
the rather large width of the slit, Ah. Particularly in the low 
density region (upstream, inside the cusped orbit) there are 
very few particles in the SPH simulation. This is a general 
resolution problem with SPH: particles follow the gas flow; 
thus low density regions near resonances, through which the 
flow moves rapidly, are often not well- resolved, whereas re- 
gions of high density such as the inner a^-ring contain many 
more particles than necessary for adequate resolution. The 
reader should note that this is even true in a symmetrized 2D 
simulation with 20000 particles. In three-dimensional simu- 
lations, especially when high density contrasts are present, 
resolution may be a much more severe problem. 

However, as the comparison shows, SPH gives sur- 
prisingly accurate solutions in the case at hand. Thus we 
may confidently move to fully self-gravitating simulations in 
which the real strengths of SPH will be used. 



5 DISCUSSION 

We have studied the structure of quasi-stationary shock 
fronts in realistic barred galaxy models with an ILR, and 
have found that both off-axis and on-axis shocks can re- 



sult in the same rotating potential, depending on the sound 
speed of the isothermal fluid. The off-axis shocks at low 
sound speeds had previously been assumed to be the rule 
in such potentials, whereas on-axis shocks had been found 
in potentials too shallow to have an ILR. We have veri- 
fied in one case that the gas flow resulting from our two- 
dimensional SPH simulation agrees approximately with that 
found in a grid-based simulation of Athanassoula (1992). 

If the quasi-equlibrium flow that results in a given 
galaxy model depends on the fluid parameters in the hydro- 
dynamical calculation — here, the sound speed — then it 
may equally depend on the entire prescription of modelling 
the ISM: whether its complicated multi-phase structure is 
better represented as particle-like or in terms of an ideal 
fluid. This offers the opportunity of learning about the 'cor- 
rect' description of interstellar fluids from detailed studies 
of galactic flows. 

The critical sound speed that we have found to divide 
the two shock regimes is around ~ 20kms -1 in our stan- 
dard model potential, which corresponds to a barred galaxy 
with a circular velocity of v c — 200kms _1 . It is somewhat 
lower for a model with a stronger bar. A dependence on 
bar strength etc. is not surprising since the occurrence of a 
shock depends on both the sound speed and the fluid flow 
velocities in the relevant bar region. 

Expressed dimensionlessly, the critical sound speed is 
~ 10% of the circular velocity. This implies that dwarf or 
Magellanic galaxies should predominantly be in the on-axis 
regime. This might have important implications for the de- 
pendence of morphology on absolute magnitude. 

For our Galaxy, the value of ~ 20kms _1 is in an inter- 
esting regime, for while the cloud velocity dispersion in the 
Galactic disk inside the solar radius is ~ 5kms~ (Clemens 
1985), the vertical cloud velocity dispersion in the Galac- 
tic bulge region is inferred to be ~ 25kms~ 1 (Bally et al. 
1988). Moreover, the effective sound speed in the ISM may 
be higher than the cloud dispersion if magnetic fields con- 
tribute significant pressure, as may be expected especially 
near the Galactic Centre (Morris 1994). 

Assuming that our Galaxy is representative, it would 
thus appear difficult to predict the type of flow pattern 
which would form in a barred galaxy. It is even unclear 
whether the flow would find a single quasi-stationary flow 
configuration. There also arises the possibility that large- 
scale star formation, by increasing the turbulent pressure in 
the disk, might itself change the structure of the flow by 
which it was initiated. Since the on-axis shock patterns are 
generally associated with larger mass inflow rates, a star- 
burst changing the pattern from off-axis to on-axis could 
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thereby prolong its own gas supply. One may even speculate 
that at early times, when star formation rates are generally 
higher, the morphology of the gas flows could have been sys- 
tematically shifted towards on-axis shock flows with higher 
mass accretion rates. If true this would have obvious rel- 
evance to HST observations of high-redshift galaxies; thus 
further work in this direction should be useful. 
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